Import dữ liệu và các thư viện cần thiết

# require pacman
if (!require("pacman")) install.packages("pacman")
## Loading required package: pacman
pacman::p_load(
  ggplot2,
  mvtnorm,
  GGally,
  corrplot,
  readxl,
  tidyverse,
  gridExtra,
  grid,
  plotly,
  ggcorrplot,
  FactoMineR,
  factoextra,
  rgl,
  scatterplot3d,
  inspectdf,
  mvnormtest,
  pastecs,
  viridis,
  psych
)
data = read_csv("Data/2017_Yellow_Taxi_Trip_Data.csv")
## New names:
## Rows: 22699 Columns: 18
## ── Column specification
## ──────────────────────────────────────────────────────── Delimiter: "," chr
## (3): tpep_pickup_datetime, tpep_dropoff_datetime, store_and_fwd_flag dbl (15):
## ...1, VendorID, passenger_count, trip_distance, RatecodeID, PULoca...
## ℹ Use `spec()` to retrieve the full column specification for this data. ℹ
## Specify the column types or set `show_col_types = FALSE` to quiet this message.
## • `` -> `...1`
data = data[, -c(1)]
head(data)
## # A tibble: 6 × 17
##   VendorID tpep_pickup_datetime tpep_dropoff_datetime passenger_count
##      <dbl> <chr>                <chr>                           <dbl>
## 1        2 3/25/2017 8:55       3/25/2017 9:09                      6
## 2        1 4/11/2017 14:53      4/11/2017 15:19                     1
## 3        1 12/15/2017 7:26      12/15/2017 7:34                     1
## 4        2 5/7/2017 13:17       5/7/2017 13:48                      1
## 5        2 4/15/2017 23:32      4/15/2017 23:49                     1
## 6        2 3/25/2017 20:34      3/25/2017 20:42                     6
## # ℹ 13 more variables: trip_distance <dbl>, RatecodeID <dbl>,
## #   store_and_fwd_flag <chr>, PULocationID <dbl>, DOLocationID <dbl>,
## #   payment_type <dbl>, fare_amount <dbl>, extra <dbl>, mta_tax <dbl>,
## #   tip_amount <dbl>, tolls_amount <dbl>, improvement_surcharge <dbl>,
## #   total_amount <dbl>
str(data)
## tibble [22,699 × 17] (S3: tbl_df/tbl/data.frame)
##  $ VendorID             : num [1:22699] 2 1 1 2 2 2 2 2 2 1 ...
##  $ tpep_pickup_datetime : chr [1:22699] "3/25/2017 8:55" "4/11/2017 14:53" "12/15/2017 7:26" "5/7/2017 13:17" ...
##  $ tpep_dropoff_datetime: chr [1:22699] "3/25/2017 9:09" "4/11/2017 15:19" "12/15/2017 7:34" "5/7/2017 13:48" ...
##  $ passenger_count      : num [1:22699] 6 1 1 1 1 6 1 1 1 1 ...
##  $ trip_distance        : num [1:22699] 3.34 1.8 1 3.7 4.37 ...
##  $ RatecodeID           : num [1:22699] 1 1 1 1 1 1 1 1 1 1 ...
##  $ store_and_fwd_flag   : chr [1:22699] "N" "N" "N" "N" ...
##  $ PULocationID         : num [1:22699] 100 186 262 188 4 161 79 237 234 239 ...
##  $ DOLocationID         : num [1:22699] 231 43 236 97 112 236 241 114 249 237 ...
##  $ payment_type         : num [1:22699] 1 1 1 1 2 1 1 1 2 1 ...
##  $ fare_amount          : num [1:22699] 13 16 6.5 20.5 16.5 9 47.5 16 9 13 ...
##  $ extra                : num [1:22699] 0 0 0 0 0.5 0.5 1 1 0 0 ...
##  $ mta_tax              : num [1:22699] 0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5 ...
##  $ tip_amount           : num [1:22699] 2.76 4 1.45 6.39 0 2.06 9.86 1.78 0 2.75 ...
##  $ tolls_amount         : num [1:22699] 0 0 0 0 0 0 0 0 0 0 ...
##  $ improvement_surcharge: num [1:22699] 0.3 0.3 0.3 0.3 0.3 0.3 0.3 0.3 0.3 0.3 ...
##  $ total_amount         : num [1:22699] 16.56 20.8 8.75 27.69 17.8 ...
duplicated_rows = data[duplicated(data),]
duplicated_rows
## # A tibble: 0 × 17
## # ℹ 17 variables: VendorID <dbl>, tpep_pickup_datetime <chr>,
## #   tpep_dropoff_datetime <chr>, passenger_count <dbl>, trip_distance <dbl>,
## #   RatecodeID <dbl>, store_and_fwd_flag <chr>, PULocationID <dbl>,
## #   DOLocationID <dbl>, payment_type <dbl>, fare_amount <dbl>, extra <dbl>,
## #   mta_tax <dbl>, tip_amount <dbl>, tolls_amount <dbl>,
## #   improvement_surcharge <dbl>, total_amount <dbl>
colSums(is.na(data))
##              VendorID  tpep_pickup_datetime tpep_dropoff_datetime 
##                     0                     0                     0 
##       passenger_count         trip_distance            RatecodeID 
##                     0                     0                     0 
##    store_and_fwd_flag          PULocationID          DOLocationID 
##                     0                     0                     0 
##          payment_type           fare_amount                 extra 
##                     0                     0                     0 
##               mta_tax            tip_amount          tolls_amount 
##                     0                     0                     0 
## improvement_surcharge          total_amount 
##                     0                     0
data_type_pickup = class(data$tpep_pickup_datetime)
data_type_dropoff = class(data$tpep_dropoff_datetime)

print(paste("Data type of tpep_pickup_datetime: ", data_type_pickup))
## [1] "Data type of tpep_pickup_datetime:  character"
print(paste("Data type of tpep_dropoff_datetime: ", data_type_dropoff))
## [1] "Data type of tpep_dropoff_datetime:  character"
data$tpep_pickup_datetime <- as.POSIXct(data$tpep_pickup_datetime, format = "%m/%d/%Y %H:%M")
data$tpep_dropoff_datetime <- as.POSIXct(data$tpep_dropoff_datetime, format = "%m/%d/%Y %H:%M")

data_type_pickup = class(data$tpep_pickup_datetime)
data_type_dropoff = class(data$tpep_dropoff_datetime)

print(paste("Data type of tpep_pickup_datetime: ", data_type_pickup))
## [1] "Data type of tpep_pickup_datetime:  POSIXct"
## [2] "Data type of tpep_pickup_datetime:  POSIXt"
print(paste("Data type of tpep_dropoff_datetime: ", data_type_dropoff))
## [1] "Data type of tpep_dropoff_datetime:  POSIXct"
## [2] "Data type of tpep_dropoff_datetime:  POSIXt"

Thông tin về dữ liệu và phương pháp nghiên cứu

Thông tin về tập dữ liệu
- ‘VendorID’ - Mã cho biết nhà cung cấp TPEP đã cung cấp hồ sơ. (1=‘Creative Mobile Technologies, LLC’, 2=‘VeriFone Inc’)
- ‘tpep_pickup_datetime’ - Ngày và giờ khi đồng hồ được bật.
- ‘tpep_dropoff_datetime’ - Ngày và giờ khi đồng hồ được ngắt.
- ‘passenger_count’ - Số lượng hành khách trên xe.
- ‘trip_distance’ - Khoảng cách chuyến đi đã trôi qua tính bằng dặm do đồng hồ taxi báo cáo.
- ‘RatecodeID’ - Mã giá cước cuối cùng có hiệu lực vào cuối chuyến đi. (1=‘Giá chuẩn’, 2=‘JFK’, 3=‘Newark’, 4=‘Nassau of Westchester’, 5=‘Giá đã thương lượng’, 6=‘Đi theo nhóm’, 99=‘NULL/không xác định’)
- ‘store_and_fwd_flag’ - Cờ này cho biết liệu hồ sơ chuyến đi có được lưu trong bộ nhớ xe trước khi gửi đến nhà cung cấp hay không, còn gọi là “lưu trữ và chuyển tiếp”, vì xe không có kết nối với máy chủ. (Y=‘lưu trữ và chuyển tiếp chuyến đi’, N=‘không phải là chuyến đi lưu trữ và chuyển tiếp’)
- ‘PULocationID’ - Khu vực taxi TLC mà đồng hồ tính cước đã được bật
- ‘DOLocationID’ - Khu vực taxi TLC mà đồng hồ tính cước đã được tắt
- ‘payment_type’ - Mã số biểu thị cách hành khách thanh toán cho chuyến đi. (0= ‘Flex Fare trip’, 1= ‘Credit card’, 2= ‘Cash’, 3= ‘No charge’, 4= ‘Dispute’, 5= ‘Unknown’, 6= ‘Voided trip’)
- ‘fare_amount’ - Giá vé theo thời gian và khoảng cách được tính theo đồng hồ. Để biết thêm thông tin về các cột sau
- ‘extra’ - Các khoản phụ phí và phụ thu khác.
- ‘mta_tax’ - Thuế được tự động kích hoạt dựa trên mức giá theo đồng hồ đang sử dụng.
- ‘tip_amount’ - Số tiền boa – Trường này được tự động điền cho tiền boa bằng thẻ tín dụng. Tiền boa bằng tiền mặt không được bao gồm.
- ‘tolls_amount’ - Tổng số tiền của tất cả các khoản phí cầu đường đã thanh toán trong chuyến đi.
- ‘improvement_surcharge’ - Phụ phí cải thiện được đánh giá cho các chuyến đi tại thời điểm hạ cờ. Phụ phí cải thiện bắt đầu được áp dụng vào năm 2015.
- ‘total_amount’ - Tổng số tiền tính cho hành khách. Không bao gồm tiền boa bằng tiền mặt.

Ghi chú về Bộ dữ liệu: Bộ dữ liệu chứa 22699 hàng và 17 cột, không có giá trị bị thiếu hoặc hàng trùng lặp. Tất cả các cột số đều là số dương.

Phân tích dữ liệu

eda = data
eda = eda[, -c(8,9)]
# Convert vendorID 1='Creative Mobile Technologies, LLC', 2='VeriFone Inc'
eda$VendorID <- ifelse(eda$VendorID == 1, 'Creative Mobile Technologies, LLC', 'VeriFone Inc')
# Convert RatecodeID 1='Standard rate', 2='JFK', 3='Newark', 4='Nassau or Westchester', 5='Negotiated fare', 6='Group ride', 99='NULL/Unknown'
eda$RatecodeID <- ifelse(eda$RatecodeID == 1, 'Standard rate', ifelse(eda$RatecodeID == 2, 'JFK', ifelse(eda$RatecodeID == 3, 'Newark', ifelse(eda$RatecodeID == 4, 'Nassau or Westchester', ifelse(eda$RatecodeID == 5, 'Negotiated fare', ifelse(eda$RatecodeID == 6, 'Group ride', 'NULL/Unknown'))))))
# Convert store_and_fwd_flag Y='store and forward trip', N='not a store and forward trip'
eda$store_and_fwd_flag <- ifelse(eda$store_and_fwd_flag == 'Y', 'store and forward trip', 'not a store and forward trip')
# Convert payment_type 0='Flex Fare trip', 1='Credit card', 2='Cash', 3='No charge', 4='Dispute', 5='Unknown', 6='Voided trip'
eda$payment_type <- ifelse(eda$payment_type == 0, 'Flex Fare trip', ifelse(eda$payment_type == 1, 'Credit card', ifelse(eda$payment_type == 2, 'Cash', ifelse(eda$payment_type == 3, 'No charge', ifelse(eda$payment_type == 4, 'Dispute', ifelse(eda$payment_type == 5, 'Unknown', 'Voided trip'))))))
# Convert to factor
eda$VendorID <- as.factor(eda$VendorID)
eda$RatecodeID <- as.factor(eda$RatecodeID)
eda$store_and_fwd_flag <- as.factor(eda$store_and_fwd_flag)
eda$payment_type <- as.factor(eda$payment_type)
str(eda)
## tibble [22,699 × 15] (S3: tbl_df/tbl/data.frame)
##  $ VendorID             : Factor w/ 2 levels "Creative Mobile Technologies, LLC",..: 2 1 1 2 2 2 2 2 2 1 ...
##  $ tpep_pickup_datetime : POSIXct[1:22699], format: "2017-03-25 08:55:00" "2017-04-11 14:53:00" ...
##  $ tpep_dropoff_datetime: POSIXct[1:22699], format: "2017-03-25 09:09:00" "2017-04-11 15:19:00" ...
##  $ passenger_count      : num [1:22699] 6 1 1 1 1 6 1 1 1 1 ...
##  $ trip_distance        : num [1:22699] 3.34 1.8 1 3.7 4.37 ...
##  $ RatecodeID           : Factor w/ 6 levels "JFK","Nassau or Westchester",..: 6 6 6 6 6 6 6 6 6 6 ...
##  $ store_and_fwd_flag   : Factor w/ 2 levels "not a store and forward trip",..: 1 1 1 1 1 1 1 1 1 1 ...
##  $ payment_type         : Factor w/ 4 levels "Cash","Credit card",..: 2 2 2 2 1 2 2 2 1 2 ...
##  $ fare_amount          : num [1:22699] 13 16 6.5 20.5 16.5 9 47.5 16 9 13 ...
##  $ extra                : num [1:22699] 0 0 0 0 0.5 0.5 1 1 0 0 ...
##  $ mta_tax              : num [1:22699] 0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5 ...
##  $ tip_amount           : num [1:22699] 2.76 4 1.45 6.39 0 2.06 9.86 1.78 0 2.75 ...
##  $ tolls_amount         : num [1:22699] 0 0 0 0 0 0 0 0 0 0 ...
##  $ improvement_surcharge: num [1:22699] 0.3 0.3 0.3 0.3 0.3 0.3 0.3 0.3 0.3 0.3 ...
##  $ total_amount         : num [1:22699] 16.56 20.8 8.75 27.69 17.8 ...
# Extract number columns
num_cols <- eda %>% select_if(is.numeric) %>% names()
num_cols
## [1] "passenger_count"       "trip_distance"         "fare_amount"          
## [4] "extra"                 "mta_tax"               "tip_amount"           
## [7] "tolls_amount"          "improvement_surcharge" "total_amount"
# stat.desc 
stats = stat.desc(eda[,num_cols])
formatted_stats <- format(stats, scientific = FALSE)
formatted_stats
##              passenger_count  trip_distance     fare_amount           extra
## nbr.val      22699.000000000 22699.00000000  22699.00000000 22699.000000000
## nbr.null        33.000000000   148.00000000      6.00000000 11921.000000000
## nbr.na           0.000000000     0.00000000      0.00000000     0.000000000
## min              0.000000000     0.00000000   -120.00000000    -1.000000000
## max              6.000000000    33.96000000    999.99000000     4.500000000
## range            6.000000000    33.96000000   1119.99000000     5.500000000
## sum          37279.000000000 66129.29000000 295691.46000000  7565.000000000
## median           1.000000000     1.61000000      9.50000000     0.000000000
## mean             1.642319045     2.91331292     13.02662937     0.333274594
## SE.mean          0.008530566     0.02424748      0.08790406     0.003073748
## CI.mean.0.95     0.016720495     0.04752673      0.17229798     0.006024756
## var              1.651819029    13.34565969    175.39798725     0.214458441
## std.dev          1.285231119     3.65317118     13.24379052     0.463096579
## coef.var         0.782570916     1.25395770      1.01667056     1.389534599
##                       mta_tax     tip_amount   tolls_amount
## nbr.val      22699.0000000000 22699.00000000 22699.00000000
## nbr.null        90.0000000000  8057.00000000 21525.00000000
## nbr.na           0.0000000000     0.00000000     0.00000000
## min             -0.5000000000     0.00000000     0.00000000
## max              0.5000000000   200.00000000    19.10000000
## range            1.0000000000   200.00000000    19.10000000
## sum          11291.5000000000 41670.40000000  7094.38000000
## median           0.5000000000     1.35000000     0.00000000
## mean             0.4974448214     1.83578131     0.31254152
## SE.mean          0.0002619441     0.01858882     0.00928710
## CI.mean.0.95     0.0005134284     0.03643536     0.01820335
## var              0.0015574852     7.84350752     1.95779403
## std.dev          0.0394649873     2.80062627     1.39921193
## coef.var         0.0793354069     1.52557729     4.47688334
##              improvement_surcharge   total_amount
## nbr.val           22699.0000000000  22699.0000000
## nbr.null              6.0000000000      4.0000000
## nbr.na                0.0000000000      0.0000000
## min                  -0.3000000000   -120.3000000
## max                   0.3000000000   1200.2900000
## range                 0.6000000000   1320.5900000
## sum                6799.5000000000 370232.0900000
## median                0.3000000000     11.8000000
## mean                  0.2995506410     16.3105022
## SE.mean               0.0001040259      0.1068439
## CI.mean.0.95          0.0002038979      0.2094213
## var                   0.0002456347    259.1229160
## std.dev               0.0156727376     16.0972953
## coef.var              0.0523208283      0.9869282

Có một số thông tin nổi bật từ bảng thống kê tóm tắt này. Rõ ràng có một số giá trị ngoại lệ trong một số biến, như tip_amount ($200) và total_amount ($1.200). Ngoài ra, một số biến, như mta_tax, dường như gần như không đổi trong toàn bộ dữ liệu, ta không mong đợi chúng mang lại khả năng dự đoán cao.

plot_all_densities <- function(data, include_cols = num_cols , fill_color = "skyblue", alpha_value = 0.3, base_size = 15) {
  plot_list <- list()
  for (var in names(data[, num_cols])) {
    p = ggplot(data, aes_string(x = var)) +
      geom_histogram(fill = "skyblue", color = "blue", alpha = alpha_value, stat='count') +  
      theme_minimal(base_size = base_size) + 
      labs(title = var,
           x = "",
           y = "Frequency",
           ) +
      geom_rug(sides = "b")
     plot_list[[var]] <- p
  }
  title_grob <- textGrob("Distribution Of Numerical Columns", gp = gpar(fontsize = 20, fontface = "bold"))
  grid.arrange(title_grob, grobs=plot_list, ncol = 3, nrow = 3)
}

plot_all_densities(eda)
## Warning: `aes_string()` was deprecated in ggplot2 3.0.0.
## ℹ Please use tidy evaluation idioms with `aes()`.
## ℹ See also `vignette("ggplot2-in-packages")` for more information.
## This warning is displayed once every 8 hours.
## Call `lifecycle::last_lifecycle_warnings()` to see where this warning was
## generated.
## Warning in geom_histogram(fill = "skyblue", color = "blue", alpha = alpha_value, : Ignoring unknown parameters: `binwidth`, `bins`, and `pad`
## Ignoring unknown parameters: `binwidth`, `bins`, and `pad`
## Ignoring unknown parameters: `binwidth`, `bins`, and `pad`
## Ignoring unknown parameters: `binwidth`, `bins`, and `pad`
## Ignoring unknown parameters: `binwidth`, `bins`, and `pad`
## Ignoring unknown parameters: `binwidth`, `bins`, and `pad`
## Ignoring unknown parameters: `binwidth`, `bins`, and `pad`
## Ignoring unknown parameters: `binwidth`, `bins`, and `pad`
## Ignoring unknown parameters: `binwidth`, `bins`, and `pad`

# countplot for factor columns
plot_all_countplots <- function(data, include_cols = c("VendorID", "RatecodeID", "store_and_fwd_flag", "payment_type"), fill_color = "skyblue", alpha_value = 0.8, base_size = 15) {
  plot_list <- list()
  for (var in include_cols) {
    p = ggplot(data, aes_string(x = var)) +
      geom_bar(fill = "aquamarine4", alpha = alpha_value) +  
      theme_minimal(base_size = base_size) + 
      theme(axis.text.x = element_text(angle = 45, hjust = 1)) +
      theme(
        plot.title = element_text(size = 16), 
        axis.title.x = element_text(size = 16), 
        axis.title.y = element_text(size = 16), 
        axis.text.x = element_text(size = 16), 
        axis.text.y = element_text(size = 16),
        legend.title = element_text(size = 16), 
        legend.text = element_text(size = 16)   
    ) +
      labs(title = var,
           x = "",
           y = "Frequency",
           ) +
      geom_text(stat='count', aes(label=..count..), vjust=-0.5)
    plot_list[[var]] <- p
  }
  title_grob <- textGrob("Countplot For Factor Columns", gp = gpar(fontsize = 20, fontface = "bold"))
  grid.arrange(title_grob, grobs=plot_list, ncol = 2, nrow = 2)
}
plot_all_countplots(eda)
## Warning: The dot-dot notation (`..count..`) was deprecated in ggplot2 3.4.0.
## ℹ Please use `after_stat(count)` instead.
## This warning is displayed once every 8 hours.
## Call `lifecycle::last_lifecycle_warnings()` to see where this warning was
## generated.

eda$duration <- as.numeric(difftime(eda$tpep_dropoff_datetime, eda$tpep_pickup_datetime, units = "mins"))
head(eda$duration)
## [1] 14 26  8 31 17  8
num_cols <- eda %>% select_if(is.numeric) %>% names()
num_cols
##  [1] "passenger_count"       "trip_distance"         "fare_amount"          
##  [4] "extra"                 "mta_tax"               "tip_amount"           
##  [7] "tolls_amount"          "improvement_surcharge" "total_amount"         
## [10] "duration"
boxplot(eda[, num_cols], rotate = 45, col = "skyblue", main = "Boxplot of numerical columns")

head(sort(unique(eda$trip_distance)), 10)
##  [1] 0.00 0.01 0.02 0.03 0.04 0.05 0.06 0.07 0.08 0.09

Khoảng cách được ghi lại với độ chính xác cao. Tuy nhiên, có thể các chuyến đi có khoảng cách bằng không nếu hành khách gọi taxi rồi đổi ý. Bên cạnh đó, có đủ giá trị bằng không trong dữ liệu để gây ra vấn đề không?

sum(eda$trip_distance == 0)
## [1] 148

148 trong số ~23.000 chuyến đi là không đáng kể. Ta có thể quy cho nó giá trị 0,01, nhưng nó không có nhiều tác động đến mô hình. Do đó, cột trip_distance sẽ không bị ảnh hưởng đối với các giá trị ngoại lệ.

summary(eda$fare_amount)
##    Min. 1st Qu.  Median    Mean 3rd Qu.    Max. 
## -120.00    6.50    9.50   13.03   14.50  999.99

Phạm vi giá trị trong cột fare_amount rất lớn và các giá trị cực trị cũng không có nhiều ý nghĩa.
* Với các giá trị thấp: Giá trị âm là có vấn đề. Giá trị bằng 0 có thể hợp lệ nếu taxi ghi lại chuyến đi đã bị hủy ngay lập tức.
* Với các giá trị cao: Số tiền cước tối đa trong tập dữ liệu này gần 1.000 đô la, điều này có vẻ rất khó xảy ra trong một chuyến taxi. Giá trị cao cho tính năng này có thể được giới hạn dựa trên trực giác và số liệu thống kê. Phạm vi tứ phân vị (IQR) là 8 đô la. Công thức chuẩn Q3 + (1,5 * IQR) cho ra kết quả là 26,50 đô la. Có vẻ không phù hợp với mức giới hạn cước tối đa. Trong trường hợp này, chúng ta sẽ sử dụng hệ số 6 * IQR, với kết quả ở mức 6*IQR là 62,5 đô la

Với các giá trị nhỏ hơn 0: Ta quy về 0

eda$fare_amount[eda$fare_amount < 0] <- 0
min(eda$fare_amount)
## [1] 0

Tiếp theo gán giá trị lớn nhất là Q3 + (6 * IQR).

# Tiếp theo gán giá trị lớn nhất là `Q3 + (6 * IQR)`.  
Q3 <- quantile(eda$fare_amount, 0.75)
IQR <- IQR(eda$fare_amount)
upper_bound <- Q3 + (6 * IQR)
eda$fare_amount[eda$fare_amount > upper_bound] <- upper_bound
summary(eda$fare_amount)
##    Min. 1st Qu.  Median    Mean 3rd Qu.    Max. 
##     0.0     6.5     9.5    12.9    14.5    62.5
summary(eda$duration)
##    Min. 1st Qu.  Median    Mean 3rd Qu.    Max. 
##  -17.00    7.00   11.00   17.01   18.00 1440.00
eda$duration[eda$duration < 0] <- 0
min(eda$duration)
## [1] 0
Q3 <- quantile(eda$duration, 0.75)
IQR <- IQR(eda$duration)
upper_bound <- Q3 + (6 * IQR)
eda$duration[eda$duration > upper_bound] <- upper_bound
summary(eda$duration)
##    Min. 1st Qu.  Median    Mean 3rd Qu.    Max. 
##    0.00    7.00   11.00   14.44   18.00   84.00
summary(eda$total_amount)
##    Min. 1st Qu.  Median    Mean 3rd Qu.    Max. 
## -120.30    8.75   11.80   16.31   17.80 1200.29

Cột total_amount cũng có các giá trị ngoại lệ ở cả hai extreme value dưới và trên.
* Giá trị thấp: Giá trị âm không hợp lý. Ta sẽ gán tất cả giá trị âm bằng 0.
* Giá trị cao: Gán các giá trị cao theo cùng phương pháp đã gán các giá trị ngoại lệ cho giá vé: Q3 + (6 * IQR).

eda$total_amount[eda$total_amount < 0] <- 0
min(eda$total_amount)
## [1] 0
Q3 <- quantile(eda$total_amount, 0.75)
IQR <- IQR(eda$total_amount)
upper_bound <- Q3 + (6 * IQR)
eda$total_amount[eda$total_amount > upper_bound] <- upper_bound
summary(eda$total_amount)
##    Min. 1st Qu.  Median    Mean 3rd Qu.    Max. 
##    0.00    8.75   11.80   16.12   17.80   72.10
boxplot(eda[, num_cols], rotate = 45, col = "skyblue", main = "Boxplot of numerical columns")

Việc sử dụng taxi thay đổi như thế nào trong một tháng?

eda$day_of_month <- as.numeric(format(eda$tpep_pickup_datetime, "%d"))
taxi_usage <- table(eda$day_of_month)
taxi_usage <- as.data.frame(taxi_usage)
colnames(taxi_usage) <- c("Day", "Count")
ggplot(taxi_usage, aes(x = Day, y = Count, group=1)) +
  geom_point(color = "coral4") +
  geom_line(color = "coral3") +
  labs(
    title = "Figure 1: How does taxi usage change over a month",
    x = "Day of the Month",
    y = "Number of Taxi Rides"
  ) +
  theme_minimal()

Chúng tôi quan sát thấy một mô hình sử dụng taxi tương đối ổn định trong suốt mỗi ngày trong tháng. Tuy nhiên, có sự sụt giảm đáng kể về số chuyến đi vào cuối tháng. Có vẻ như họ đã sử dụng hết ngân sách đi lại của mình vào đầu và giữa tháng.

Tháng nào có số chuyến đi taxi nhiều nhất trong năm?

eda$month <- as.numeric(format(eda$tpep_pickup_datetime, "%m"))

taxi_usage_month <- table(eda$month)
taxi_usage_month <- as.data.frame(taxi_usage_month)
colnames(taxi_usage_month) <- c("Month", "Count")
ggplot(taxi_usage_month, aes(x = factor(Month), y = Count, fill = factor(Month))) +
  geom_bar(stat = "identity") +
  geom_text(aes(label = Count), vjust = -0.5, color = "black") +
  scale_fill_viridis_d() +
  labs(
    title = "Figure 2: Which month has the most taxi rides in the year",
    x = "Month",
    y = "Number of Taxi Rides"
  ) +
  theme_minimal() +
  theme(legend.position = "none")

Nhà cung cấp nào có số chuyến đi taxi nhiều nhất?

taxi_usage_vendor <- table(eda$VendorID)
taxi_usage_vendor <- as.data.frame(taxi_usage_vendor)
colnames(taxi_usage_vendor) <- c("Vendor", "Count")
ggplot(taxi_usage_vendor, aes(x = Vendor, y = Count, fill = Vendor)) +
  geom_bar(stat = "identity") +
  geom_text(aes(label = Count), vjust = -0.5, color = "black") +
  scale_fill_viridis_d() +
  labs(
    title = "Figure 3: Which vendor has the most taxi rides",
    x = "Vendor",
    y = "Number of Taxi Rides"
  ) +
  theme_minimal() +
  theme(legend.position = "none")

Nhận xét: Có sự khác biệt về khoảng cách chuyến đi giữa các nhà cung cấp, nhưng không quá lớn.

Ảnh hưởng của factor đến total_amount

grid.arrange(
  ggplot(eda, aes(x = VendorID, y = trip_distance, fill = VendorID)) +
    geom_boxplot() +
    labs(
      title = "Figure 4: Is there a difference in trip distance between vendors",
      x = "Vendor",
      y = "Trip Distance"
    ) +
    theme_minimal(),
  ggplot(eda, aes(x = RatecodeID, y = total_amount, fill = RatecodeID)) +
    geom_boxplot() +
    labs(
      title = "Figure 5: Influence of RatecodeID on total amount",
      x = "RatecodeID",
      y = "Total Amount"
    ) +
    theme_minimal(),
  ggplot(eda, aes(x = store_and_fwd_flag, y = total_amount, fill = store_and_fwd_flag)) +
    geom_boxplot() +
    labs(
      title = "Figure 6: Influence of store_and_fwd_flag on total amount",
      x = "store_and_fwd_flag",
      y = "Total Amount"
    ) +
    theme_minimal(),
  ggplot(eda, aes(x = payment_type, y = total_amount, fill = payment_type)) +
    geom_boxplot() +
    labs(
      title = "Figure 7: Influence of payment_type on total amount",
      x = "payment_type",
      y = "Total Amount"
    ) +
    theme_minimal(),
  ncol = 2
)

Doanh thu cao nhất theo ngày trong tuần, tháng, quý

eda$day_of_week <- as.numeric(format(eda$tpep_pickup_datetime, "%u"))
eda$quarter <- as.numeric(format(eda$tpep_pickup_datetime, "%m")) %/% 4 + 1
daily_revenue <- eda %>%
  group_by(day_of_week) %>%
  summarise(total_amount = sum(fare_amount, na.rm = TRUE)) %>%
  arrange(day_of_week)

monthly_revenue <- eda %>%
  group_by(month) %>%
  summarise(total_amount = sum(fare_amount, na.rm = TRUE)) %>%
  arrange(month)

quarterly_revenue <- eda %>%
  group_by(quarter) %>%
  summarise(total_amount = sum(fare_amount, na.rm = TRUE)) %>%
  arrange(quarter)

highest_daily_revenue <- daily_revenue %>%
  filter(total_amount == max(total_amount, na.rm = TRUE))

highest_monthly_revenue <- monthly_revenue %>%
  filter(total_amount == max(total_amount, na.rm = TRUE))

highest_quarterly_revenue <- quarterly_revenue %>%
  filter(total_amount == max(total_amount, na.rm = TRUE))




print(paste("Ngày có doanh thu cao nhất trong tuần là thứ 5", 
            ".Với số tiền là:", highest_daily_revenue$total_amount))
## [1] "Ngày có doanh thu cao nhất trong tuần là thứ 5 .Với số tiền là: 45180.6"
print(paste("Tháng có doanh thu cao nhất trong năm: ", highest_monthly_revenue$month, 
            "Với số tiền là: ", highest_monthly_revenue$total_amount))
## [1] "Tháng có doanh thu cao nhất trong năm:  5 Với số tiền là:  26823.01"
print(paste("Qúy có doanh thu cao nhất trong năm là quý", highest_quarterly_revenue$quarter, 
            "Với số tiền là", highest_quarterly_revenue$total_amount))
## [1] "Qúy có doanh thu cao nhất trong năm là quý 2 Với số tiền là 99782.81"

Heatmap correlation

correlation_matrix <- cor(eda[, num_cols])
ggcorrplot(correlation_matrix, lab = TRUE, method = "square", title = "Figure 8: Correlation heatmap of numerical columns")

Kiểm định MANOVA

manova_results <- manova(cbind(trip_distance, fare_amount, duration, total_amount) ~ VendorID, data = eda)
summary(manova_results)
##              Df     Pillai approx F num Df den Df  Pr(>F)  
## VendorID      1 0.00037978   2.1555      4  22694 0.07131 .
## Residuals 22697                                            
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Kiểm định ANOVA

print("Kiểm định ANOVA cho total_amount")
## [1] "Kiểm định ANOVA cho total_amount"
anova_results <- aov(total_amount ~ VendorID, data = eda)
summary(anova_results)
##                Df  Sum Sq Mean Sq F value Pr(>F)
## VendorID        1      17   16.62     0.1  0.752
## Residuals   22697 3780099  166.55
anova_results <- aov(duration ~ VendorID, data = eda)
print("Kiểm định ANOVA cho duration")
## [1] "Kiểm định ANOVA cho duration"
summary(anova_results)
##                Df  Sum Sq Mean Sq F value Pr(>F)
## VendorID        1      53   53.27   0.379  0.538
## Residuals   22697 3190306  140.56
anova_results <- aov(fare_amount ~ VendorID, data = eda)
print("Kiểm định ANOVA cho fare_amount")
## [1] "Kiểm định ANOVA cho fare_amount"
summary(anova_results)
##                Df  Sum Sq Mean Sq F value Pr(>F)
## VendorID        1       3    2.76   0.025  0.875
## Residuals   22697 2522099  111.12
anova_results <- aov(trip_distance ~ VendorID, data = eda)
print("Kiểm định ANOVA cho trip_distance")
## [1] "Kiểm định ANOVA cho trip_distance"
summary(anova_results)
##                Df Sum Sq Mean Sq F value Pr(>F)
## VendorID        1     18   18.34   1.375  0.241
## Residuals   22697 302901   13.35

Phân tích dữ liệu bằng phương pháp PCA

data.pca <- PCA(eda[, num_cols], graph=F)
eig.val <- get_eigenvalue(data.pca)
eig.val
##         eigenvalue variance.percent cumulative.variance.percent
## Dim.1  4.496876686      44.96876686                    44.96877
## Dim.2  1.643460437      16.43460437                    61.40337
## Dim.3  1.002979679      10.02979679                    71.43317
## Dim.4  0.974999128       9.74999128                    81.18316
## Dim.5  0.674774045       6.74774045                    87.93090
## Dim.6  0.588084482       5.88084482                    93.81174
## Dim.7  0.316920815       3.16920815                    96.98095
## Dim.8  0.219676290       2.19676290                    99.17772
## Dim.9  0.077175198       0.77175198                    99.94947
## Dim.10 0.005053241       0.05053241                   100.00000

Scree Plot

fviz_eig(data.pca, addlabels = TRUE, ylim = c(0, 100)) +
  labs(title = "Scree Plot", x = "Principal Component", y = "Percentage of Variance") +
  theme_minimal()

Loading

loadings <- as.data.frame(data.pca$var$coord)
loadings$Variable <- rownames(loadings)
loadings
##                               Dim.1        Dim.2        Dim.3       Dim.4
## passenger_count        0.0157836723 -0.007700298  0.963301835  0.26357361
## trip_distance          0.9315639649  0.061576239  0.015392126 -0.01897864
## fare_amount            0.9664719333  0.039491509  0.016399632 -0.03472738
## extra                  0.1393237374  0.156149640 -0.266129046  0.93972532
## mta_tax               -0.1762988861  0.893668718  0.018737831 -0.05074638
## tip_amount             0.6606644141 -0.016653396 -0.042364596 -0.05371105
## tolls_amount           0.7269844610 -0.105736408  0.006839232 -0.04751297
## improvement_surcharge  0.0002759055  0.889661305  0.031990056 -0.11326393
## total_amount           0.9854914266  0.043918282 -0.000493820 -0.01014541
## duration               0.8414739681  0.100696906  0.021965234 -0.01513325
##                             Dim.5              Variable
## passenger_count        0.04359601       passenger_count
## trip_distance         -0.15577176         trip_distance
## fare_amount           -0.15424734           fare_amount
## extra                  0.04069800                 extra
## mta_tax               -0.10255222               mta_tax
## tip_amount             0.58928165            tip_amount
## tolls_amount           0.30264521          tolls_amount
## improvement_surcharge  0.20360820 improvement_surcharge
## total_amount          -0.03550355          total_amount
## duration              -0.36204872              duration
p1 <- ggplot(loadings, aes(x = reorder(Variable, Dim.1), y = Dim.1)) +
  geom_bar(stat = "identity", fill = "cornsilk4") +
  coord_flip() +
  labs(title = "Loadings on PC1", x = "Variables", y = "Loadings") +
  theme_minimal() +
  theme(
        plot.title = element_text(size = 14), 
        axis.title.x = element_text(size = 14), 
        axis.title.y = element_text(size = 14), 
        axis.text.x = element_text(size = 14), 
        axis.text.y = element_text(size = 14),
        legend.title = element_text(size = 14), 
        legend.text = element_text(size = 14)   
    )

p2 <- ggplot(loadings, aes(x = reorder(Variable, Dim.2), y = Dim.2)) +
  geom_bar(stat = "identity", fill = "coral3") +
  coord_flip() +
  labs(title = "Loadings on PC2", x = "Variables", y = "Loadings") +
  theme_minimal() +
  theme(
        plot.title = element_text(size = 14), 
        axis.title.x = element_text(size = 14), 
        axis.title.y = element_text(size = 14), 
        axis.text.x = element_text(size = 14), 
        axis.text.y = element_text(size = 14),
        legend.title = element_text(size = 14), 
        legend.text = element_text(size = 14)   
    )

p3 <- ggplot(loadings, aes(x = reorder(Variable, Dim.3), y = Dim.3)) +
  geom_bar(stat = "identity", fill = "aquamarine4") +
  coord_flip() +
  labs(title = "Loadings on PC3", x = "Variables", y = "Loadings") +
  theme_minimal() + 
  theme(
        plot.title = element_text(size = 14), 
        axis.title.x = element_text(size = 14), 
        axis.title.y = element_text(size = 14), 
        axis.text.x = element_text(size = 14), 
        axis.text.y = element_text(size = 14),
        legend.title = element_text(size = 14), 
        legend.text = element_text(size = 14)   
    )

grid.arrange(p1, p2, p3, nrow = 2, ncol = 2)

Contribution of Variables to Principal Components

contrib_PC1 <- as.data.frame(data.pca$var$contrib[,1])
colnames(contrib_PC1) <- c("Contribution")
contrib_PC1$Variable <- rownames(contrib_PC1)

contrib_PC2 <- as.data.frame(data.pca$var$contrib[,2])
colnames(contrib_PC2) <- c("Contribution")
contrib_PC2$Variable <- rownames(contrib_PC2)

contrib_PC3 <- as.data.frame(data.pca$var$contrib[,3])
colnames(contrib_PC3) <- c("Contribution")
contrib_PC3$Variable <- rownames(contrib_PC3)

p1 = ggplot(contrib_PC1, aes(x = reorder(Variable, Contribution), y = Contribution)) +
  geom_bar(stat = "identity", fill = "skyblue") +
  coord_flip() +
  labs(title = "Contribution of Variables to Dim 1", x = "", y = "Contribution (%)") +
  theme_minimal()

p2 = ggplot(contrib_PC2, aes(x = reorder(Variable, Contribution), y = Contribution)) +
  geom_bar(stat = "identity", fill = "salmon") +
  coord_flip() +
  labs(title = "Contribution of Variables to Dim 2", x = "", y = "Contribution (%)") +
  theme_minimal()

p3 = ggplot(contrib_PC3, aes(x = reorder(Variable, Contribution), y = Contribution)) +
  geom_bar(stat = "identity", fill = "purple") +
  coord_flip() +
  labs(title = "Contribution of Variables to Dim 3", x = "", y = "Contribution (%)") +
  theme_minimal()

grid.arrange(p1,
             p2,
             p3,
             nrow=1)

Biplot

p1 <-  fviz_pca_var(data.pca,
               col.var = "contrib",
               gradient.cols = c("#00AFBB", "#E7B800", "#FC4E07"),
               repel = TRUE)

p2 <- fviz_pca_var(data.pca, axes = c(2, 3),
             col.var = "contrib",
             gradient.cols = c("#00AFBB", "#E7B800", "#FC4E07"),
             repel = TRUE)
grid.arrange(p1, p2, nrow = 1)

fviz_pca_biplot(data.pca, axes = c(1, 2), geom = "point", col.var = "contrib")

fviz_pca_biplot(data.pca, axes = c(2, 3), geom = "point", col.var = "contrib")

# merge biplot pc1 pc2 pc3
biplot_pc1_pc2 <- fviz_pca_biplot(data.pca, axes = c(1, 2), geom = "point", col.var = "contrib")
biplot_pc2_pc3 <- fviz_pca_biplot(data.pca, axes = c(2, 3), geom = "point", col.var = "contrib")
grid.arrange(biplot_pc1_pc2, biplot_pc2_pc3, nrow = 1)

# pca 3d
explained_variance <- data.pca$eig[, 2] 
pca_3d <- as.data.frame(data.pca$ind$coord)
pca_3d$VendorID <- data$VendorID

plot_ly(pca_3d, x = ~Dim.1, y = ~Dim.2, z = ~Dim.3, color = ~VendorID, colors = c("#00AFBB", "salmon")) %>%
  add_markers() %>%
  layout(scene = list(xaxis = list(title = paste0('Dim1 (', round(explained_variance[1], 2), '%)')),
                      yaxis = list(title = paste0('Dim2 (', round(explained_variance[2], 2), '%)')),
                      zaxis = list(title = paste0('Dim3 (', round(explained_variance[3], 2), '%)'))),
         title = "3D PCA Plot of Psychological Data")
eda
## # A tibble: 22,699 × 20
##    VendorID           tpep_pickup_datetime tpep_dropoff_datetime passenger_count
##    <fct>              <dttm>               <dttm>                          <dbl>
##  1 VeriFone Inc       2017-03-25 08:55:00  2017-03-25 09:09:00                 6
##  2 Creative Mobile T… 2017-04-11 14:53:00  2017-04-11 15:19:00                 1
##  3 Creative Mobile T… 2017-12-15 07:26:00  2017-12-15 07:34:00                 1
##  4 VeriFone Inc       2017-05-07 13:17:00  2017-05-07 13:48:00                 1
##  5 VeriFone Inc       2017-04-15 23:32:00  2017-04-15 23:49:00                 1
##  6 VeriFone Inc       2017-03-25 20:34:00  2017-03-25 20:42:00                 6
##  7 VeriFone Inc       2017-05-03 19:04:00  2017-05-03 20:03:00                 1
##  8 VeriFone Inc       2017-08-15 17:41:00  2017-08-15 18:03:00                 1
##  9 VeriFone Inc       2017-02-04 16:17:00  2017-02-04 16:29:00                 1
## 10 Creative Mobile T… 2017-11-10 15:20:00  2017-11-10 15:40:00                 1
## # ℹ 22,689 more rows
## # ℹ 16 more variables: trip_distance <dbl>, RatecodeID <fct>,
## #   store_and_fwd_flag <fct>, payment_type <fct>, fare_amount <dbl>,
## #   extra <dbl>, mta_tax <dbl>, tip_amount <dbl>, tolls_amount <dbl>,
## #   improvement_surcharge <dbl>, total_amount <dbl>, duration <dbl>,
## #   day_of_month <dbl>, month <dbl>, day_of_week <dbl>, quarter <dbl>

Phân tích nhân tố (FA)

Thông thường, trước khi tiến hành phân tích thành phần chính (PCA) và EFA, ta phải kiểm tra xem dữ liệu ta thu thập được có phù hợp cho phân tích hay không. Hai phương pháp kiểm định thường dùng để làm điều này là kiểm định KMO và Bartlett
## Kiểm định KMO

kmo_fa <- KMO(eda[, num_cols])
kmo_fa
## Kaiser-Meyer-Olkin factor adequacy
## Call: KMO(r = eda[, num_cols])
## Overall MSA =  0.64
## MSA for each item = 
##       passenger_count         trip_distance           fare_amount 
##                  0.70                  0.95                  0.61 
##                 extra               mta_tax            tip_amount 
##                  0.17                  0.33                  0.50 
##          tolls_amount improvement_surcharge          total_amount 
##                  0.60                  0.35                  0.62 
##              duration 
##                  0.96

Một bộ dữ liệu được xem là phù hợp để phân tích EFA thường có giá trị KMO Test tối thiểu là 0.5. Trong trường hợp này, giá trị KMO Test là 0.64, vì vậy ta có thể tiến hành phân tích nhân tố cho bộ dữ liệu.

Kiểm định Bartlett

cortest.bartlett(eda[, num_cols])
## R was not square, finding R from data
## $chisq
## [1] 214695.6
## 
## $p.value
## [1] 0
## 
## $df
## [1] 45

Bartlett Test có giá trị 19334.49, bậc tự do là 45 và trị số p (p value =0) nhỏ hơn <5%, có nghĩa là quan hệ giữa các biến (items) đủ lớn để sử dụng phân tích EFA.

Xác định số lượng các nhân tố chính rút ra

Trong phân tích EFA, căn cứ để xác định các nhân tố chính được rút ra là sử dụng giá trị của Eigenvalue. Theo tiêu chuẩn của Kaiser thì nhân tố chính được rút ra phải có Eigenvalue > 1. Một tiêu chuẩn khác ít nghiêm ngặt hơn đó là tiêu chuẩn của Jolliffe, theo Jolliffe thì các nhân tố có Eigenvalue > 0.7 có thể được chọn. Trong bộ dữ liệu này, số lượng nhân tố chính được rút ra dựa vào tiêu chuẩn của Kaiser.

fviz_eig(data.pca, addlabels = TRUE, ylim = c(0, 100), n=10) +
  labs(title = "Scree Plot", x = "Principal Component", y = "Percentage of Variance") +
  theme_minimal()

Xác định các biến cấu thành nhân tố được rút ra, đặt tên nhân tố

Ở trên, sử dụng Eigenvalue theo tiêu chuẩn Kaiser ta đã trích ra được 3 nhân tố chính, để biết các nhân tố này được cấu thành từ những biến (items) nào ta sử dụng phép xoay Varimax.
print.psych(pc2, cut = 3, sort = TRUE) # Loại bỏ các biến có hệ số tải (Factor Loading) bé hơn 0.4

pc2 = principal(eda[, num_cols], nfactors = 3, rotate = "varimax")
print.psych(pc2, cut = 0.4, sort = TRUE)
## Principal Components Analysis
## Call: principal(r = eda[, num_cols], nfactors = 3, rotate = "varimax")
## Standardized loadings (pattern matrix) based upon correlation matrix
##                       item   RC1   RC2   RC3   h2    u2 com
## total_amount             9  0.98             0.97 0.027 1.0
## fare_amount              3  0.97             0.94 0.064 1.0
## trip_distance            2  0.93             0.87 0.128 1.0
## duration                10  0.84             0.72 0.281 1.0
## tolls_amount             7  0.72             0.54 0.460 1.1
## tip_amount               6  0.66             0.44 0.561 1.0
## mta_tax                  5        0.90       0.83 0.170 1.1
## improvement_surcharge    8        0.89       0.79 0.207 1.0
## passenger_count          1              0.96 0.93 0.072 1.0
## extra                    4                   0.11 0.885 1.8
## 
##                        RC1  RC2  RC3
## SS loadings           4.48 1.64 1.02
## Proportion Var        0.45 0.16 0.10
## Cumulative Var        0.45 0.61 0.71
## Proportion Explained  0.63 0.23 0.14
## Cumulative Proportion 0.63 0.86 1.00
## 
## Mean item complexity =  1.1
## Test of the hypothesis that 3 components are sufficient.
## 
## The root mean square of the residuals (RMSR) is  0.07 
##  with the empirical chi square  9489.04  with prob <  0 
## 
## Fit based upon off diagonal values = 0.97
  • Theo kết quả trên ta thấy với nhân tố thứ nhất, các biến cấu thành nên nhân tố này là trip_distance, fare_amount, duration, total_amount, tolls_amount’, ‘tip_amount’. Các biến này đều liên quan đến các khía cạnh chi phí và khoảng cách của chuyến đi, bao gồm chi phí cơ bản, chi phí phụ (như phí cầu đường và tiền tip), và các yếu tố thời gian và khoảng cách của chuyến đi. Vì vậy, ta có thể đặt tên cho nhân tố thứ nhất này là “Chi phí và Khoảng cách Chuyến đi”.
  • Với nhân tố thứ hai, các biến cấu thành nên nhân tố này là ‘mta_tax’, ‘improvement_surcharge’. Các biến này đều liên quan đến các khoản phí và thuế bổ sung được tính thêm vào chi phí chuyến đi. Vì vậy, ta có thể đặt tên cho nhân tố thứ hai này là “Phí và Thuế Bổ sung.
  • Với nhân tố thứ ba, các biến cấu thành nên nhân tố này là ‘passenger_count’. Biến này liên quan đến số lượng hành khách trong mỗi chuyến đi. Vì vậy, ta có thể đặt tên cho nhân tố thứ ba này là “Số lượng Hành khách”.

Kiểm định Cronbach Alpha cho thang đo

Để kiểm tra thang đo nhân tố được rút ra được tạo thành từ các biến có phù hợp hay không ta sử dụng kiểm định Cronbach Alpha. Chẳng hạn để kiểm định xem nhân tố thứ nhất đặt tên là “Chi phí và khoảng cách” được cấu thành từ 6 biến có phù hợp không, ta sử dụng lệnh sau

psych::alpha(eda[, c('trip_distance', 'fare_amount', 'duration', 'total_amount', 'tolls_amount', 'tip_amount')])
## Number of categories should be increased  in order to count frequencies.
## 
## Reliability analysis   
## Call: psych::alpha(x = eda[, c("trip_distance", "fare_amount", "duration", 
##     "total_amount", "tolls_amount", "tip_amount")])
## 
##   raw_alpha std.alpha G6(smc) average_r S/N     ase mean  sd median_r
##       0.86      0.93    0.95      0.67  12 0.00046  8.1 6.6     0.63
## 
##     95% confidence boundaries 
##          lower alpha upper
## Feldt     0.86  0.86  0.87
## Duhachek  0.86  0.86  0.86
## 
##  Reliability if an item is dropped:
##               raw_alpha std.alpha G6(smc) average_r  S/N alpha se var.r med.r
## trip_distance      0.84      0.90    0.94      0.64  8.9  0.00055 0.037  0.62
## fare_amount        0.77      0.89    0.91      0.62  8.2  0.00080 0.031  0.63
## duration           0.82      0.92    0.95      0.68 10.8  0.00043 0.037  0.63
## total_amount       0.78      0.89    0.90      0.61  7.8  0.00093 0.032  0.58
## tolls_amount       0.89      0.93    0.96      0.73 13.8  0.00047 0.039  0.79
## tip_amount         0.88      0.94    0.95      0.76 15.8  0.00044 0.029  0.79
## 
##  Item statistics 
##                   n raw.r std.r r.cor r.drop  mean   sd
## trip_distance 22699  0.92  0.92  0.92   0.91  2.91  3.7
## fare_amount   22699  0.98  0.96  0.98   0.96 12.90 10.5
## duration      22699  0.90  0.84  0.80   0.81 14.44 11.9
## total_amount  22699  0.98  0.98  1.00   0.97 16.12 12.9
## tolls_amount  22699  0.64  0.74  0.67   0.62  0.31  1.4
## tip_amount    22699  0.61  0.69  0.64   0.56  1.84  2.8

Một thang đo được xem là hợp lý nếu giá trị Cronbach Alpha > 0.7. Trong bộ dữ liệu này, giá trị Cronbach Alpha là ~0,86 nên có thể nói rằng thang đo “Chi phí và khoảng cách” được cấu hình thành từ 7 biến (items) là hợp lý.Tương tự, thực hiện kiểm tra thang đo cho 2 nhân tố còn lại.

psych::alpha(eda[, c('mta_tax', 'improvement_surcharge')])
## 
## Reliability analysis   
## Call: psych::alpha(x = eda[, c("mta_tax", "improvement_surcharge")])
## 
##   raw_alpha std.alpha G6(smc) average_r S/N    ase mean    sd median_r
##        0.6      0.77    0.63      0.63 3.4 0.0031  0.4 0.025     0.63
## 
##     95% confidence boundaries 
##          lower alpha upper
## Feldt     0.59   0.6  0.61
## Duhachek  0.60   0.6  0.61
## 
##  Reliability if an item is dropped:
##                       raw_alpha std.alpha G6(smc) average_r S/N alpha se var.r
## mta_tax                    1.58      0.63     0.4      0.63 1.7       NA     0
## improvement_surcharge      0.25      0.63     0.4      0.63 1.7       NA     0
##                       med.r
## mta_tax                0.63
## improvement_surcharge  0.63
## 
##  Item statistics 
##                           n raw.r std.r r.cor r.drop mean    sd
## mta_tax               22699  0.97   0.9  0.72   0.63  0.5 0.039
## improvement_surcharge 22699  0.80   0.9  0.72   0.63  0.3 0.016
## 
## Non missing response frequency for each item
##                       -0.5 -0.3 0 0.3 0.5 miss
## mta_tax                  0    0 0   0   1    0
## improvement_surcharge    0    0 0   1   0    0
psych::alpha(eda[, c('passenger_count', 'extra', 'mta_tax', 'improvement_surcharge')])
## Number of categories should be increased  in order to count frequencies.
## Warning in psych::alpha(eda[, c("passenger_count", "extra", "mta_tax", "improvement_surcharge")]): Some items were negatively correlated with the first principal component and probably 
## should be reversed.  
## To do this, run the function again with the 'check.keys=TRUE' option
## Some items ( passenger_count ) were negatively correlated with the first principal component and 
## probably should be reversed.  
## To do this, run the function again with the 'check.keys=TRUE' option
## 
## Reliability analysis   
## Call: psych::alpha(x = eda[, c("passenger_count", "extra", "mta_tax", 
##     "improvement_surcharge")])
## 
##   raw_alpha std.alpha G6(smc) average_r  S/N    ase mean   sd median_r
##    -0.0036      0.35    0.41      0.12 0.53 0.0057 0.69 0.34     0.02
## 
##     95% confidence boundaries 
##          lower alpha upper
## Feldt    -0.03     0  0.02
## Duhachek -0.01     0  0.01
## 
##  Reliability if an item is dropped:
##                       raw_alpha std.alpha G6(smc) average_r   S/N alpha se
## passenger_count         2.2e-02     0.486   0.503     0.240 0.947  0.00211
## extra                   7.4e-05     0.439   0.479     0.207 0.783  0.00079
## mta_tax                -5.7e-03     0.032   0.023     0.011 0.033  0.00641
## improvement_surcharge  -5.1e-03     0.037   0.027     0.013 0.039  0.00644
##                         var.r   med.r
## passenger_count       0.11324  0.0508
## extra                 0.13330 -0.0014
## mta_tax               0.00067 -0.0014
## improvement_surcharge 0.00109 -0.0063
## 
##  Item statistics 
##                           n raw.r std.r   r.cor  r.drop mean    sd
## passenger_count       22699 0.939  0.42 -0.0094 -0.0068 1.64 1.285
## extra                 22699 0.335  0.47  0.0592 -0.0043 0.33 0.463
## mta_tax               22699 0.047  0.72  0.7196  0.0185 0.50 0.039
## improvement_surcharge 22699 0.042  0.72  0.7157  0.0306 0.30 0.016